function [Zp,Z,miu,niu,miu0,niu0] = ...
    updateAllZ(C,gamma,eta,alpha,Zp,Z,miu_old,niu_old,miu0_old,niu0_old)
% Computing Multi-view Optimal Transport
% Update miu_p, niu_p, miu_0 and niu_0 by Newton's method. 
% X:d*n*v, C:n*m*v, A:d*m*v, alpha:1*1*v, Zp:n*m*v, Z:n*m, miu:n*1*v, niu:1*m*v, miu0:n*1, niu0:1*m

%% Precomputation
[numsample,numanchor,numview] = size(C);
% alpha = reshape(alpha,1,1,numview);
b = alpha.^2 * eta + gamma;
b0 = sum(alpha.^2) * eta;

%% Check input
if nargin < 5 || isempty(Zp)
    Zp = repmat(1/numsample/numanchor,numsample,numanchor,numview);%
    Z = repmat(1/numsample/numanchor,numsample,numanchor);%
end
if nargin < 7 || isempty(miu_old)
    D = C - pagemtimes(repmat(Z,1,1,numview) , alpha.^2) * eta * 2;
%     niu_old = zeros(1,numanchor,numview);
%     miu_old = - mean(pagemtimes(Zp , b)*2 + D + repmat(niu_old,numsample,1),2);
%     niu_old = - mean(pagemtimes(Zp , b)*2 + D + repmat(miu_old,1,numanchor),1);
    miu_old = - mean(pagemtimes(Zp , b) + D/2,2);
    niu_old = - mean(pagemtimes(Zp , b) + D/2,1);
    Zp = pagemtimes(max(-(D + repmat(miu_old,1,numanchor) + repmat(niu_old,numsample,1)),0) , 0.5 ./ b);
    
    D0 = - sum(pagemtimes(Zp , alpha.^2),3) * eta * 2;
%     niu0_old = zeros(1,numanchor);
%     miu0_old = - mean(Z * b0*2 + D0 + repmat(niu0_old,numsample,1),2);
%     niu0_old = - mean(Z * b0*2 + D0 + repmat(miu0_old,1,numanchor),1);
    miu0_old = - mean(Z * b0 + D0/2,2);
    niu0_old = - mean(Z * b0 + D0/2,1);
    Z = max(-(D0 + repmat(miu0_old,1,numanchor) + repmat(niu0_old,numsample,1)),0) / b0 / 2;
end

%% Preset
e = 1e-3;

iter = 0;
max_iter = 20;
max_iter0 = 20;

miu = miu_old;
niu = niu_old;
miu0 = miu0_old;
niu0 = niu0_old;

%% Iteriation...
while (iter <= max_iter) && ~(iter >= 1 && sum(abs(sum(Z,2) - 1),'all') <= e)
    iter = iter + 1;
    
    % Update miu_p, niu_p by Newton's method. 
    D = C - pagemtimes(repmat(Z,1,1,numview) , alpha.^2) * eta * 2;
    Zp = pagemtimes(max(-(D + repmat(miu,1,numanchor) + repmat(niu,numsample,1)),0) , 0.5 ./ b);
    iter0 = 0;
    while (iter0 <= max_iter0) && ~(iter0 >= 1 && sum(abs(sum(Zp,2) - 1),'all') <= e)
        iter0 = iter0 + 1;
        miu = miu + pagemtimes((sum(Zp,2) - 1/numsample) ./ max(sum(sign(Zp),2),1e-3) , b) * 2;
        Zp = pagemtimes(max(-(D + repmat(miu,1,numanchor) + repmat(niu,numsample,1)),0) , 0.5 ./ b);
    end
    
    iter0 = 0;
    while (iter0 <= max_iter0) && ~(iter0 >= 1 && sum(abs(sum(Zp,1) - 1),'all') <= e)
        iter0 = iter0 + 1;
        niu = niu + pagemtimes((sum(Zp,1) - 1/numanchor) ./ max(sum(sign(Zp),1),1e-3) , b) * 2;
        Zp = pagemtimes(max(-(D + repmat(miu,1,numanchor) + repmat(niu,numsample,1)),0) , 0.5 ./ b);
    end
    
    % Update miu_0, niu_0 by Newton's method. 
    D0 = - sum(pagemtimes(Zp , alpha.^2),3) * eta * 2;
    Z = max(-(D0 + repmat(miu0,1,numanchor) + repmat(niu0,numsample,1)),0) / b0 / 2;
    iter0 = 0;
    while (iter0 <= max_iter0) && ~(iter0 >= 1 && sum(abs(sum(Z,2) - 1),'all') <= e)
        iter0 = iter0 + 1;        
        miu0 = miu0 + (sum(Z,2) - 1/numsample) ./ max(sum(sign(Z),2),1e-3) * b0 * 2;
        Z = max(-(D0 + repmat(miu0,1,numanchor) + repmat(niu0,numsample,1)),0) / b0 / 2;
    end
    
    iter0 = 0;
    while (iter0 <= max_iter0) && ~(iter0 >= 1 && sum(abs(sum(Z,1) - 1),'all') <= e)
        iter0 = iter0 + 1;        
        niu0 = niu0 + (sum(Z,1) - 1/numanchor) ./ max(sum(sign(Z),1),1e-3) * b0 * 2;
        Z = max(-(D0 + repmat(miu0,1,numanchor) + repmat(niu0,numsample,1)),0) / b0 / 2;
    end
end




